module duty
  !
  ! Dicom Format Utilities
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  use suty
  implicit none
  private
  integer,parameter::METAHDRSIZMAX=4096
  INTEGER, PARAMETER :: SEEK_SET = 0, SEEK_CUR = 1, SEEK_END = 2
  integer,parameter::UID_MAXSTR = 64
  character*(*),parameter::STR_UID_CREATOR="777.777.0.0.0"
  integer,parameter::UID_CREATOR            = 1
  integer,parameter::UID_SOP_INSTANCE       = 2
  integer,parameter::UID_MEDIA_INSTANCE     = 3
  integer,parameter::UID_FRAME              = 4
  integer,parameter::UID_STUDY              = 5
  integer,parameter::UID_SERIES             = 6
  public write_dcm
  public image_info,study_info,patient_info,etc_info
  type study_info
     integer year,month,day,hour,minute,second
     character*32 id
     character*128 description
     character*32 injectionTime
     character*32 part
     character*32 rphar
     real dose
     real halflife
     real positronFrac
     character*128 ofile
     character*128 input
     integer duration
  end type study_info
  type patient_info
     character*32 name
     character*32 id
     character*32 dob
     character*32 sex
     integer age
     real weight,height
  end type patient_info
  type etc_info
     character*128 institution
     character*128 operator
     character*128 derivDesc
  end type etc_info
  type image_info
     real slice_width
     real slice_separation
     integer number_of_slices
     integer*2 width,height
     real dx,dy
     real bias,intcpt
     integer hdrsz
     type(study_info)   sinf
     type(patient_info) pinf
     type(etc_info)     einf
  end type image_info
  integer p    ! output to pointer
  integer unit ! output to file
  character*(4)::MG_DICM="DICM"
  integer,parameter::DCM_HDR_MAX=8192

  !
  type(image_info),pointer::ifo
  type(study_info),pointer::sinf
  type(patient_info),pointer::pinf
  type(etc_info),pointer::einf

  character*128 imgdate,imgtime

  !
contains

  !
  integer function hashSDBM(str)
    character*(*),intent(in)::str
    integer i
    hashSDBM=0
    do i=1,len_trim(str)
       hashSDBM=ichar(str(i:i))+ishft(hashSDBM,6)+ishft(hashSDBM,16)-hashSDBM
    end do
    hashSDBM=iand(hashSDBM,int(z'FFFFFFFF',8))
  end function hashSDBM

  !
  integer function mkUID(type,uid)
    integer,intent(in)::type
    character*(*),intent(out)::uid
    integer study_uid,series_uid,instance_uid
    integer utc
    character*128 str
    
    call system_clock(count=utc)
    if(utc==-huge(utc)) utc=777

    write(str,100) trim(pinf%name),trim(pinf%id),&
         sinf%year,sinf%month,sinf%day,sinf%hour,sinf%minute,sinf%second
    str=adjustl(str)
100 format(2a,i4,5(i2.2))
    
    study_uid=hashSDBM(str)
    if(study_uid==182208422) then
       str=trim(i2a(utc))
       study_uid=hashSDBM(str)
    end if

    series_uid=hashSDBM(sinf%input)
    instance_uid=hashSDBM(sinf%ofile)

    select case(type)
       case(UID_CREATOR)
          uid=STR_UID_CREATOR
       case(UID_SOP_INSTANCE,UID_MEDIA_INSTANCE)
          uid=STR_UID_CREATOR//"."//trim(i2a(utc))//"."//trim(i2a(study_uid))//"."&
               //trim(i2a(series_uid))//"."//trim(i2a(instance_uid))
       case(UID_FRAME,UID_STUDY)
          uid=STR_UID_CREATOR//"."//trim(i2a(utc))//"."//trim(i2a(study_uid))

       case(UID_SERIES)
          uid=STR_UID_CREATOR//"."//trim(i2a(utc))//"."//trim(i2a(study_uid))//"."//trim(i2a(series_uid))

       case default
          uid=STR_UID_CREATOR//"."//trim(i2a(utc))//"."//trim(i2a(study_uid))
       end select

       mkUid=len_trim(uid)
       if(mkUid>UID_MAXSTR) then
          write(*,*) "Inappropriate UID length."
       end if
     end function mkUID

  !
  integer function write_metahdr() result(istat)
    integer*1,allocatable::metahdr(:)
    integer p0,unit0
    integer*1::hdr(128)=0_1
    integer*2::OB(7)=[int(z'0002',2),&
         int(z'0001',2),&
         int(z'424f',2),&
         int(z'0000',2),&
         int(z'0002',2),&
         int(z'0000',2),&
         int(z'0100',2)]
    character*128 buf
!    integer istat
    integer rewrf,begin,endp
    integer len
    integer*4 dw

    allocate(metahdr(METAHDRSIZMAX))
    p0=p
    unit0=unit
    unit=0
    p=loc(metahdr)
    p=cpmemf1(p,loc(hdr),sizeof(hdr))
    p=cpmemf1(p,loc(MG_DICM),sizeof(MG_DICM))

    rewrf=p
    istat=write_element(int8(z'0002'),int8(z'0000'),"??",4,loc(hdr))
    begin=p
    
    p=cpmemf1(p,loc(OB),sizeof(OB))
    
    buf="1.2.840.10008.5.1.4.1.1.20" ! NM modality
    istat=write_element(int8(z'0002'),int8(z'0002'),"UI",len_trim(buf),loc(buf))

    len=mkUid(UID_MEDIA_INSTANCE,buf)
    istat=write_element(int8(z'0002'),int8(z'0003'),"UI",len,loc(buf))

    buf="1.2.840.10008.1.2.1"
    istat=write_element(int8(z'0002'),int8(z'0010'),"UI",len_trim(buf),loc(buf))

    buf="0.0.0.0"
    istat=write_element(int8(z'0002'),int8(z'0012'),"UI",len_trim(buf),loc(buf))

    endp=p
    dw=int(endp-begin,4)
    istat=write_element_at(int8(z'0002'),int8(z'0000'),"UL",4,loc(dw),rewrf)

    if(unit0>0) then
       len=p-loc(metahdr)
       write(unit0,iostat=istat) metahdr(1:len)
    elseif(p0>0) then
       call cpmems1(p0,loc(metahdr),p-1)
    end if
    deallocate(metahdr)
    p=p0
    unit=unit0
  end function write_metahdr

  !
  integer function write_element_at(group,element,vr,length,pdata,pos) result(istat)
    integer,intent(in)::group,element,length,pdata
    character*2,intent(in)::vr
    integer,intent(in)::pos
    integer pcur
    pcur=p
    p=pos
    istat=write_element(group,element,vr,length,pdata)
    p=pcur
  end function write_element_at

  !
  integer function write_element(group,element,vr,length,pdata) result(istat)
    integer,intent(in)::group,element,length,pdata
    character*2,intent(in)::vr
    integer len
    logical make_even
    integer*1 b
    integer*1 pad
    integer*2 w
    pointer(pp,b)
    istat=0
    if(iand(length,1)/=0) then
       len=length+1
       make_even=.true.
    else
       len=length
       make_even=.false.
    end if
    if(unit>0) then
       write(unit,iostat=istat) int(group,kind=2),int(element,kind=2),vr,int(len,kind=2)
    elseif(p>0) then
       w=int(group,kind=2)
       p=cpmemf2(p,loc(w),1)
       w=int(element,kind=2)
       p=cpmemf2(p,loc(w),1)
       p=cpmemf2(p,loc(vr),1)
       w=int(len,kind=2)
       p=cpmemf2(p,loc(w),1)
    end if
    if(istat/=0) return
    if(length>0) then
       do pp=pdata,pdata+length-1
          if(unit>0) then
             write(unit,iostat=istat) b
             if(istat/=0) return
          elseif(p>0) then
             p=cpmemf1(p,loc(b),1)
          end if
       end do
       if(make_even) then
          select case(vr)
          case("UI")
             pad=0_1
          case default
             pad=32_1
          end select
          if(unit>0) then
             write(unit,iostat=istat) pad
             if(istat/=0) return
          elseif(p>0) then
             p=cpmemf1(p,loc(pad),1)
          end if
       end if
    end if
  end function write_element

  !
  integer function write_dcm(ounit,ifo0,pimg,len,pimgout) result(istat)
    integer,intent(in)::ounit
    type(image_info),intent(in),target::ifo0
    integer,intent(in)::pimg,len
    integer,intent(in),optional::pimgout ! output
    integer*1,allocatable::buf(:)
    integer*2 img(*)
    pointer(pimg,img)

    ifo=>ifo0
    pinf=>ifo0%pinf
    sinf=>ifo0%sinf
    einf=>ifo0%einf
    imgdate=""
    imgtime=""
    call date_and_time(date=imgdate,time=imgtime)
    if(present(pimgout)) then
       if(pimgout<=0) then
          allocate(buf(DCM_HDR_MAX))
          p=loc(buf)
          istat=write_hdr()
          istat=p-loc(buf)+1
          deallocate(buf)
          return
       else
          p=pimg
       end if
       unit=-1
    else
       p=-1
       unit=ounit
    end if

    !
    istat=write_hdr()
    if(istat/=0) then
       istat=-istat
       return
    end if

    !
    if(unit>0) then
       write(unit,iostat=istat) img(1:len)
    elseif(p>0) then
       p=cpmemf1(p,loc(img),len*2) 
       istat=p-pimg+1
    end if
  end function write_dcm

  !
  integer function write_hdr() result(istat)
    use suty
    character*128 buf
    character*1,parameter::delm=achar(z'5c')
    integer len
    integer*2 w
    integer*2::w7fe0=int2(z'7fe0')
    integer*2::w0010=int2(z'0010')
    integer*2::w0000=int2(z'0000')
    character*(2)::cOW="OW"
    integer*4 dw

    !
    istat=write_metahdr()
    
    buf="DERIVED"//delm//"PRIMARY"//delm//"RECON TOMO"//delm//"EMISSION"
!    buf="DERIVED"//delm//"PRIMARY"//delm//"STATIC"//delm//"EMISSION"
    istat=write_element(int8(z'0008'),int8(z'0008'),"CS",len_trim(buf),loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0012'),"DA",len_trim(imgdate),loc(imgdate))
    istat=write_element(int8(z'0008'),int8(z'0013'),"TM",len_trim(imgtime),loc(imgtime))    
    len=mkUid(UID_CREATOR,buf)
    istat=write_element(int8(z'0008'),int8(z'0014'),"UI",len,loc(buf))
    buf="1.2.840.10008.5.1.4.1.1.20" ! NM modality
    istat=write_element(int8(z'0008'),int8(z'0016'),"UI",len_trim(buf),loc(buf))
    len=mkUid(UID_SOP_INSTANCE,buf)
    istat=write_element(int8(z'0008'),int8(z'0018'),"UI",len,loc(buf))
    write(buf,100) sinf%year,sinf%month,sinf%day
    buf=adjustl(buf)
100 format(i4,2(i2.2))
    len=len_trim(buf)
    istat=write_element(int8(z'0008'),int8(z'0020'),"DA",len,loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0021'),"DA",len,loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0022'),"DA",len,loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0023'),"DA",len_trim(imgdate),loc(imgdate))
    write(buf,200) sinf%hour,sinf%minute,sinf%second
    buf=adjustl(buf)
200 format(3(i2.2))
    len=len_trim(buf)
    istat=write_element(int8(z'0008'),int8(z'0030'),"TM",len,loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0031'),"TM",len,loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0032'),"TM",len,loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0033'),"TM",len_trim(imgtime),loc(imgtime))
    buf="1"
    istat=write_element(int8(z'0008'),int8(z'0050'),"SH",len_trim(buf),loc(buf))
    buf="NM"
    istat=write_element(int8(z'0008'),int8(z'0060'),"CS",len_trim(buf),loc(buf))
    buf="ELSEIF LABORATORY"
    istat=write_element(int8(z'0008'),int8(z'0070'),"LO",len_trim(buf),loc(buf))
    istat=write_element(int8(z'0008'),int8(z'0080'),"LO",len_trim(einf%institution),loc(einf%institution))

    buf="JESUSCHRIST" ! referring physicians name
    istat=write_element(int8(z'0008'),int8(z'0090'),"PN",len_trim(buf),loc(buf))

    istat=write_element(int8(z'0008'),int8(z'1030'),"LO",len_trim(sinf%description),loc(sinf%description))
    istat=write_element(int8(z'0008'),int8(z'103e'),"LO",len_trim(sinf%description),loc(sinf%description))
    istat=write_element(int8(z'0008'),int8(z'1070'),"PN",len_trim(einf%operator),loc(einf%operator))

    istat=write_element(int8(z'0008'),int8(z'2111'),"ST",len_trim(einf%derivDesc),loc(einf%derivDesc))    

    istat=write_element(int8(z'0010'),int8(z'0010'),"PN",len_trim(pinf%name),loc(pinf%name))
    istat=write_element(int8(z'0010'),int8(z'0020'),"LO",len_trim(pinf%id),loc(pinf%id))
    buf=trim(strip(pinf%dob))
    istat=write_element(int8(z'0010'),int8(z'0030'),"DA",len_trim(buf),loc(buf))

    buf=""
    istat=write_element(int8(z'0010'),int8(z'0032'),"TM",len_trim(buf),loc(buf)) !<<

    buf=pinf%sex(1:1)
    istat=write_element(int8(z'0010'),int8(z'0040'),"CS",len_trim(buf),loc(buf))
    buf=f2a(pinf%height,-1)
    istat=write_element(int8(z'0010'),int8(z'1020'),"DS",len_trim(buf),loc(buf))    
    buf=f2a(pinf%weight,-1)
    istat=write_element(int8(z'0010'),int8(z'1030'),"DS",len_trim(buf),loc(buf))    

    istat=write_element(int8(z'0018'),int8(z'0015'),"CS",len_trim(sinf%part),loc(sinf%part))

    istat=write_element(int8(z'0018'),int8(z'0031'),"LO",len_trim(sinf%rphar),loc(sinf%rphar))

    buf=f2a(ifo%slice_width,-1)
    istat=write_element(int8(z'0018'),int8(z'0050'),"DS",len_trim(buf),loc(buf))

    buf=f2a(ifo%slice_separation,-1)
    istat=write_element(int8(z'0018'),int8(z'0088'),"DS",len_trim(buf),loc(buf))

    buf=strip(sinf%injectionTime)
    if(len_trim(buf)==4) buf=trim(buf)//"00"
    istat=write_element(int8(z'0018'),int8(z'1072'),"DS",len_trim(buf),loc(buf))
    buf=f2a(sinf%dose,-1)
    istat=write_element(int8(z'0018'),int8(z'1074'),"DS",len_trim(buf),loc(buf))
    buf=f2a(sinf%halflife,-1)
    istat=write_element(int8(z'0018'),int8(z'1075'),"DS",len_trim(buf),loc(buf))
    buf=f2a(sinf%positronFrac,-1)
    istat=write_element(int8(z'0018'),int8(z'1076'),"DS",len_trim(buf),loc(buf)) ! positron fraction

    buf=i2a(sinf%duration)
    istat=write_element(int8(z'0018'),int8(z'1242'),"IS",len_trim(buf),loc(buf))
    
    !
    len=mkUid(UID_STUDY,buf)
    istat=write_element(int8(z'0020'),int8(z'000d'),"UI",len,loc(buf))
    len=mkUid(UID_SERIES,buf)
    istat=write_element(int8(z'0020'),int8(z'000e'),"UI",len,loc(buf))
    istat=write_element(int8(z'0020'),int8(z'0010'),"SH",len_trim(sinf%id),loc(sinf%id)) !<<
    buf="0"
    istat=write_element(int8(z'0020'),int8(z'0011'),"IS",len_trim(buf),loc(buf))
    istat=write_element(int8(z'0020'),int8(z'0012'),"IS",len_trim(buf),loc(buf))
    istat=write_element(int8(z'0020'),int8(z'0013'),"IS",len_trim(buf),loc(buf))
    len=mkUid(UID_STUDY,buf)
    istat=write_element(int8(z'0020'),int8(z'0052'),"UI",len,loc(buf))
    buf=i2a(ifo%number_of_slices)
    istat=write_element(int8(z'0020'),int8(z'1002'),"IS",len_trim(buf),loc(buf))
    istat=write_element(int8(z'0020'),int8(z'1040'),"LO",0,0)

    w=1_2
    istat=write_element(int8(z'0028'),int8(z'0002'),"US",2,loc(w))
    buf="MONOCHROME2"
    istat=write_element(int8(z'0028'),int8(z'0004'),"CS",len_trim(buf),loc(buf))
    buf=i2a(ifo%number_of_slices)
    istat=write_element(int8(z'0028'),int8(z'0008'),"IS",len_trim(buf),loc(buf))    

!< (0028,0009) AT[1] FrameIncrementPointer: (0054,0080) (4 bytes)
    
    dw=int(z'00800054',4)
    istat=write_element(int8(z'0028'),int8(z'0009'),"AT",4,loc(dw)) !<<

    istat=write_element(int8(z'0028'),int8(z'0010'),"US",2,loc(ifo%height)) ! rows
    istat=write_element(int8(z'0028'),int8(z'0011'),"US",2,loc(ifo%width))! columns
    buf=trim(adjustl(f2a(ifo%dx,-1)))//delm//trim(adjustl(f2a(ifo%dy,-1)))
    istat=write_element(int8(z'0028'),int8(z'0030'),"DS",len_trim(buf),loc(buf))
    w=int(z'0010',2)
    istat=write_element(int8(z'0028'),int8(z'0100'),"US",2,loc(w))
    w=int(z'0010',2)
    istat=write_element(int8(z'0028'),int8(z'0101'),"US",2,loc(w))
    w=int(z'000f',2)
    istat=write_element(int8(z'0028'),int8(z'0102'),"US",2,loc(w))
    w=int(z'0001',2)
    istat=write_element(int8(z'0028'),int8(z'0103'),"US",2,loc(w))

    buf=f2a(ifo%intcpt,-1)
    istat=write_element(int8(z'0028'),int8(z'1052'),"DS",len_trim(buf),loc(buf))
    buf=f2a(ifo%bias,-1)
    istat=write_element(int8(z'0028'),int8(z'1053'),"DS",len_trim(buf),loc(buf))

    w=int(ifo%number_of_slices,2)
    istat=write_element(int8(z'0054'),int8(z'0081'),"US",2,loc(w)) !<<

    len=ifo%number_of_slices*ifo%width*ifo%height*2
    if(unit>0) then
       write(unit,iostat=istat) w7fe0,w0010,cOW,w0000,int(len,kind=4)
       if(istat/=0) return
    elseif(p>0) then
       p=cpmemf1(p,loc(w7fe0),sizeof(w7fe0))
       p=cpmemf1(p,loc(w0010),sizeof(w0010))
       p=cpmemf1(p,loc(cOW),sizeof(cOW))
       p=cpmemf1(p,loc(w0000),sizeof(w0000))
       p=cpmemf1(p,loc(len),4)
    end if
  end function write_hdr

  !
  character*32 function strip(str)
    character*(*),intent(in)::str
    integer i,k
    strip=""
    k=0
    do i=1,len_trim(str)
       select case(str(i:i))
       case("0","1","2","3","4","5","6","7","8","9",".")
          k=k+1
          strip(k:k)=str(i:i)
       end select
    end do
  end function strip

  !  
end module duty
